Vamos a buscar data de repositorios online, sobretodo usaremos Open Street Map (OSM)
# Cargamos las libraries
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(ggmap)
## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
## Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
## OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
OSM tiene informacion de todo el mundo y vamos a querer bajar solo un área de interés. Vamos a ponerle como string la query de la búsqueda del lugar.
La wiki de OSM tiene la info necesaria para entender el catálogo: https://wiki.openstreetmap.org/wiki/ES:Objetos_del_mapa
# Obtenemos la bounding box del área de interes
bbox_aoi <- getbb("Comuna 10, Ciudad Autonoma de Buenos Aires, Buenos Aires, Argentina")
bbox_aoi
## min max
## x -58.53145 -58.47131
## y -34.64581 -34.60912
Chequeemos que las coordenadas de la bounding box estan ok haciendo un mapa y viendo que sea la zona del mundo que nos interesa, no sea cosa que estemos en otro hemisferio porejemplo. Podría pasar, se imaginan?
register_stadiamaps("cebdf63b-1fbd-40ac-af86-fc993f352022", write = TRUE)
## ℹ Replacing old key (cebdf63b) with new key in /Users/danielarisaro/.Renviron
mapa_aoi <- get_stadiamap(bbox = bbox_aoi,
maptype = "stamen_toner_lite",
zoom = 14)
## ℹ © Stadia Maps © Stamen Design © OpenMapTiles © OpenStreetMap contributors.
ggmap(mapa_aoi)
Ademas del area de la comuna tambien podemos usar los poligonos que las delimitan
poligono_aoi <- getbb("Comuna 10, Ciudad Autonoma de Buenos Aires, Buenos Aires, Argentina",
format_out = "sf_polygon")
Podemos hacer la superposicion del area y del poligono para ver como lucen juntos
ggmap(mapa_aoi) +
geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 1.25, inherit.aes = FALSE) +
labs(title = "Comuna 10",
subtitle = "Ciudad Autonoma de Buenos Aires",
caption = "Fuente: OpenStreetMap.") +
theme_void()
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
La informacion disponible de OSM la podemos consultar con dos
funciones: available_features() y
available_tags() Las features nos dan keys y cada key tiene
tags
all_features <- available_features()
print(all_features[1:10])
## [1] "4wd_only" "abandoned"
## [3] "abutters" "access"
## [5] "addr" "addr:city"
## [7] "addr:conscriptionnumber" "addr:country"
## [9] "addr:county" "addr:district"
available_tags('building')
## # A tibble: 101 × 2
## Key Value
## <chr> <chr>
## 1 building allotment_house
## 2 building annexe
## 3 building apartments
## 4 building bakehouse
## 5 building barn
## 6 building barracks
## 7 building beach_hut
## 8 building boathouse
## 9 building bridge
## 10 building bungalow
## # ℹ 91 more rows
Como opciones podemos descargar lineas, puntos o poligonos (y multipolígonos). Vamos de a uno?
Vamos a descargar algunos puntos en particular de las churchs
churchs_aoi <- opq(bbox_aoi)
churchs_aoi <- add_osm_feature(churchs_aoi,
key = "building",
value = c("building", "church"))
churchs_aoi <- osmdata_sf(churchs_aoi)
churchs_aoi
## Object of class 'osmdata' with:
## $bbox : -34.6458134,-58.5314494,-34.609115,-58.4713084
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 126 points
## $osm_lines : NULL
## $osm_polygons : 'sf' Simple Features Collection with 12 polygons
## $osm_multilines : NULL
## $osm_multipolygons : NULL
Me quedo con los puntos. Me podria quedar con los poligonos tambien
haciendo osm_polygons
churchs_aoi <- churchs_aoi$osm_points
A ver cuantos puntos tenemos:
dim(churchs_aoi)
## [1] 126 3
126 puntos!
ggmap(mapa_aoi) +
geom_sf(data = churchs_aoi, inherit.aes = FALSE) +
geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
labs(title = "Edificios religiosos - iglesias",
subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
caption = "Fuente: OpenStreetMap") +
theme_void() +
theme(title = element_text(size = 8, face = "bold"), #tamaño de titulo del mapa
plot.caption = element_text(face = "italic", colour = "gray35", size = 7)) #tamaño de nota al pie
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Vamos a ver exactamente cuantas churchs entran en la comuna 10. Para
eso vamos a usar la st_intersection()
churchs_aoi <- st_intersection(churchs_aoi, poligono_aoi)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
dim(churchs_aoi)
## [1] 68 3
son sesentayocho. che me parece UN MONTON. Y por qué hay puntos tan contiguos en el mapa? Pues no lo se, pero habría que ver qué entiende OSM por church
Tambien podemos descargar lineas del area of interest.
calles_aoi <- opq(bbox_aoi)
calles_aoi <- add_osm_feature(calles_aoi,
key = "highway",
value = c("motorway", "trunk", "primary", "secondary", "tertiary"))
calles_aoi <- osmdata_sf(calles_aoi)
calles_aoi
## Object of class 'osmdata' with:
## $bbox : -34.6458134,-58.5314494,-34.609115,-58.4713084
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 2701 points
## $osm_lines : 'sf' Simple Features Collection with 718 linestrings
## $osm_polygons : 'sf' Simple Features Collection with 0 polygons
## $osm_multilines : NULL
## $osm_multipolygons : NULL
Nos quedamos con las lines
calles_aoi <- calles_aoi$osm_lines
ggmap(mapa_aoi) +
geom_sf(data = calles_aoi, color = "#2b2d42", inherit.aes = FALSE, lwd = 0.7) +
geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
labs(title = "Calles principales",
subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
caption = "Fuente: OpenStreetMap") +
theme_void() +
theme(title = element_text(size = 8, face = "bold"), #tamaño de titulo del mapa
plot.caption = element_text(face = "italic", colour = "gray35", size = 7)) #tamaño de nota al pie
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Para finalizar vamos a hacer una query de poligonos (overpass query, remember). Acá me anoto una tarea que es ver bien el tema de los tags porque no los entiendo del todo: en los values estan concatenadas diferentes keys o hay una adentro de la otra?
parques_aoi <- opq(bbox_aoi)
parques_aoi <- add_osm_feature(parques_aoi,
key = "leisure",
value = c("park", "garden", "nature_reserve"))
parques_aoi <- osmdata_sf(parques_aoi)
parques_aoi
## Object of class 'osmdata' with:
## $bbox : -34.6458134,-58.5314494,-34.609115,-58.4713084
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 1879 points
## $osm_lines : 'sf' Simple Features Collection with 48 linestrings
## $osm_polygons : 'sf' Simple Features Collection with 110 polygons
## $osm_multilines : NULL
## $osm_multipolygons : 'sf' Simple Features Collection with 8 multipolygons
Ahora vamos a tomar los poligonos y multipoligonos de este objeto
parques_polygons <- parques_aoi$osm_polygons
parques_multipolygons <- parques_aoi$osm_multipolygons
dim(parques_polygons)
## [1] 110 29
ggmap(mapa_aoi) +
geom_sf(data = parques_polygons, fill = "#588157", color = NA, inherit.aes = FALSE) +
geom_sf(data = parques_multipolygons, fill = "#588157", color = NA, inherit.aes = FALSE) +
geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
labs(title = "Espacios verdes",
subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
caption = "Fuente: OpenStreetMap") +
theme_void()
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Tenemos dos elementos que son los poligonos y los multipoligonos y queremos que sea una sola! Para eso vamos a seleccionar by the id and then we will bind them together oh yeah
parques_polygons <- parques_polygons %>%
select(osm_id, name)
parques_multipolygons <- parques_multipolygons %>%
select(osm_id, name)
parques_aoi <- rbind(parques_polygons, parques_multipolygons)
Ahora intersectamos el poligono de la comuna de interes con los poligonos de los parques verdes
parques_aoi <- st_intersection(parques_aoi, poligono_aoi)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
ggmap(mapa_aoi) +
geom_sf(data = parques_aoi, fill = "#588157", color =NA, inherit.aes = FALSE) +
geom_sf(data = poligono_aoi, fill = NA, color = "#0077b6", lwd = 0.75, inherit.aes = FALSE) +
labs(title = "Espacios verdes",
subtitle = "Comuna 10, Ciudad Autonoma de Buenos Aires",
caption = "Fuente: OpenStreetMap") +
theme_void() +
theme(title = element_text(size = 8, face = "bold"), #tamaño de titulo del mapa
plot.caption = element_text(face = "italic", colour = "gray35", size = 7)) #tamaño de nota al pie
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
jajajaja me quedan los parques on top de los names pero ya fue
Vamos a calcular la superficie de las áreas verdes
parques_aoi <- parques_aoi %>%
mutate(area_m2 = st_area(parques_aoi))
summary(parques_aoi$area_m2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.55 239.33 876.10 2604.68 2479.92 18665.13
El area promedio de los espacios verdes en la comuna 10 es de 2600 m2 creo yo que es muy poco poquisimo, terrible mala performance
summarise(parques_aoi, area_m2 = sum(area_m2))
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -58.52973 ymin: -34.64361 xmax: -58.47425 ymax: -34.614
## Geodetic CRS: WGS 84
## area_m2 geometry
## 1 210979.2 [m^2] MULTIPOLYGON (((-58.4916 -3...
El área total de los parques es de 210979 m2, 0.2 km2